#define LOCAL

#include "task.h"
#include "method.h"
#include "math.h"
#include "parallel.h"
#include <stdio.h>
#include <math.h>
#include <iostream>
#include <algorithm>

#include <mkl.h>
#include <mkl_lapacke.h>

void PrintPGrad(Task& task)
{
	Method* method = task.CurrentMethod;
	Model* model = task.CurrentModel;
	Grid& grid = task.CurrentGrid;
	std::vector<Vector> P_gr(grid.cells.size());
	for(int i=0; i<grid.cells.size(); i++)		
	{
		std::vector<Vector> gradCell(5);
		gradCell = task.CurrentMethod->grad[i];
		P_gr[i] = gradCell[4];	// grad(ro*E)
		Vector Velocity(method->values[i][1], method->values[i][2], method->values[i][3]);	//ro*velocity in cell[i]
		Velocity /= method->values[i][0];		//velocity in cell[i]
		P_gr[i] += Velocity.x*(gradCell[1]-Velocity.x*gradCell[0]);			//grad(ro*E) + ro*u*grad_u				
		P_gr[i] += Velocity.y*(gradCell[2]-Velocity.y*gradCell[0]);			//grad(ro*E) + ro*u*grad_u + ro*v*grad_v
		P_gr[i] += Velocity.z*(gradCell[3]-Velocity.z*gradCell[0]);			//grad(ro*E) + ro*u*grad_u + ro*v*grad_v + ro*w*grad_w
		P_gr[i] += (-0.5)*(Velocity*Velocity)*gradCell[0];	 // grad(ro*E) - 0.5*V^2*grad_ro + ro*u*grad_u + ro*v*grad_v + ro*w*grad_w
		P_gr[i] *= model->gamma-1;					//gradient of P for cell[i]
	};/*
	for(int i=0; i<grid.cells.size(); i++)
	{
		std::vector<Vector> gradCell(5);
		gradCell = task.CurrentMethod->grad[i];
		P_gr[i] = gradCell[4];	// grad(ro*E)
		Vector Velocity(method->values[i][1], method->values[i][2], method->values[i][3]);	//ro*velocity in cell[i]
		Velocity /= method->values[i][0];		//velocity in cell[i]
		P_gr[i] += (-1)*Velocity.x*(gradCell[1]-Velocity.x*gradCell[0]);			//grad(ro*E) - ro*u*grad_u				
		P_gr[i] += (-1)*Velocity.y*(gradCell[2]-Velocity.y*gradCell[0]);			//grad(ro*E) - ro*u*grad_u - ro*v*grad_v
		P_gr[i] += (-1)*Velocity.z*(gradCell[3]-Velocity.z*gradCell[0]);			//grad(ro*E) - ro*u*grad_u - ro*v*grad_v - ro*w*grad_w
		P_gr[i] += (-0.5)*(Velocity*Velocity)*gradCell[0];				// grad(ro*E) - 0.5*V^2*grad_ro - ro*u*grad_u - ro*v*grad_v - ro*w*grad_w
		P_gr[i] *= model->gamma-1;					//gradient of P for cell[i]
	};*/
	std::ofstream ofs("Grad_P.dat");
	//TO DO unify		
	ofs<<"VARIABLES= \"X\", \"Y\", \"Gr_P\", \n";
	ofs<<"ZONE T=\"D\"\n";
	ofs<<"N=" << grid.nodes.size() << ", E=" << grid.cells.size() <<", F=FEBLOCK, ET=QUADRILATERAL\n";
	ofs<<"VARLOCATION = (NODAL, NODAL, CELLCENTERED)\n";
	//X
	for (int i = 0; i<grid.nodes.size(); i++) {
		ofs<<grid.nodes[i].P.x<<"\n";
	};
	//Y
	for (int i = 0; i<grid.nodes.size(); i++) {
		ofs<<grid.nodes[i].P.y<<"\n";
	};
	//Gr_P
	for (int i = 0; i<grid.cells.size(); i++) {
		ofs<<P_gr[i].x<<"\n";
	};
	for (int i = 0; i<grid.cells.size(); i++) {
		ofs<<grid.cells[i].Nodes[0]<<" "
			<<grid.cells[i].Nodes[1]<<" "
			<<grid.cells[i].Nodes[2]<<" "
			<<grid.cells[i].Nodes[3]<<"\n";
	};
	ofs.close();

};
int main(int argc, char *argv[])
{	
	/*
	int dims[2]={4,3};
	int periods[2]= {0,0}; 
	int reorder = 0;
	int coords[2];
	int cartrank;
	int x, y;
	struct {
		int nx,ny;
	} g;

	g.nx = 13;
	g.ny = 7;

	int dx = g.nx / dims[0];
	int dy = g.ny / dims[1];
	//Map cell coords to process coords		
	const int dim_max = 1000;
	int xdispl[dim_max];
	int ydispl[dim_max];
	int r = g.nx % dims[0];
	xdispl[0] = 0;
	for (int i = 1; i<=dims[0]; i++) {
		xdispl[i] = xdispl[i-1] + dx;
		if (r > 0) {
			xdispl[i]++;
			r--;
		};
	}
	r = g.ny % dims[1];
	ydispl[0] = 0;
	for (int i = 1; i<=dims[1]; i++) {
		ydispl[i] = ydispl[i-1] + dy;
		if (r > 0) {
			ydispl[i]++;
			r--;
		};
	}

	std::vector<int> c_coords(2);
	c_coords[0] = 10;
	c_coords[1] = 5;
	//TO DO modify distribution like in decomposition
	for (int i = 1; i<=dims[0]; i++) {
		if (c_coords[0] < xdispl[i]) {
			coords[0] = i-1;
			break;
		};
	};
	for (int i = 1; i<=dims[1]; i++) {
		if (c_coords[1] < ydispl[i]) {
			coords[1] = i-1;
			break;
		};
	};
	std::cout<<coords[0]<<" "<<coords[1]<<"\n";
	return 0;
	*/
	//U_average calculation
	/*Face face;
	double L = 1.0;
	double Re_l_cr = 3.2e5;
	CellBoundarySolution CBS(100, face);
	CBS.U_inf = 70.0;
	double my = 1.5e-5;
	CBS.ny_kinetic = my;
	double Re = L*CBS.U_inf/CBS.ny_kinetic;
	double h_l = 5*L/sqrt(Re);	//laminar h
	double h_t = 0.37*L/exp(0.2 * log(Re));	//turbulent h
	//CBS.h = h_l;	
	CBS.h = h_t;
	double c = 1.0;
	double alpha = 0.74*c*c*exp(0.8*log(Re_l_cr));			
	double gamma0 = 0.108;
	CBS.sigma = 2*CBS.ny_kinetic*CBS.U_inf/(CBS.h*CBS.h);
	CBS.h = h_t;
	CBS.InitAproach();
	CBS.Solve(c, alpha, gamma0, 0.01, 500);
	double I = CBS.norm()/CBS.ny_kinetic;
	CBS.sigma = 2*CBS.ny_kinetic*(CBS.U_inf-I)/(CBS.h*CBS.h);
	CBS.ComputeU();	//old type of U
	for(int i=0; i<CBS.size ; i++) CBS.U[i] /= CBS.U[CBS.size-1];
	CBS.UndimY();
	CBS.ToTecPlot("hi_4.dat", "U_turb_4.dat");
	return 0;*/
	//////////Cf test////////////
	Face face;
	double L = 1.0;
	double Re_l_cr = 3.2e5;
	CellBoundarySolution CBS(100, face);
	CBS.KARMAN = 4.1;
	CBS.U_inf = 30;
	double my = 1.51e-5;
	CBS.ny_kinetic = my;
	double c = 1.0;			//optimal value for continue U model
	double alpha = 9.8*c*c*exp(0.5*log(Re_l_cr));
	double gamma0 = 0.0051;		//optimal value for continue U model
	std::ofstream ofs("Cf.dat"); 
	std::ofstream ofs_Sliht("Cf_theory.dat");
	double norm;
	for(int i=99; i<100; i++)
	{
		std::cout << i << "\n";
		//CBS.U_inf = CBS.U_inf + i;
		//CBS.U_inf++;
		double Re = L*CBS.U_inf/CBS.ny_kinetic;
		CBS.h = CBS.height(L, Re_l_cr);
		CBS.sigma = 2*CBS.ny_kinetic*CBS.U_inf/(CBS.h*CBS.h);
		double Re_din = CBS.h*sqrt(CBS.h*CBS.sigma)/CBS.ny_kinetic;
		double Re_cr = sqrt(alpha)/c;
		if(Re_din>Re_cr)
		{
			CBS.InitAproach();
			CBS.Solve(c, alpha, gamma0, 0.01, 500);
			//CBS.TrivialSolution();
			//for (int j=0; j<CBS.hi.size(); j++) CBS.hi[j] *= 10;
			CBS.ComputeU_2();
			CBS.UndimU_plus();
			CBS.UndimYplus();
		}else{
			CBS.TrivialSolution();
			CBS.sigma = 0;
			CBS.ComputeU_1();
		};
		for(int j=1; j<CBS.y.size(); j++) CBS.y[j] = log10(CBS.y[j]);
		CBS.ToTecPlot();
		for(int j=0; j<CBS.U.size(); j++) CBS.U[j] = CBS.U[j]*(CBS.U_inf - CBS.U[j]);	// U(U_inf - u) function
		double W_specific = CBS.ComputeU_av_h()*CBS.h;
		double Cf = 2*W_specific/(CBS.U_inf*CBS.U_inf*L);			//for flat plate (Slihting)
		ofs << log10(Re) << '\t' << 10*log10(1000*Cf) << '\n';
		
		if(Re<=Re_l_cr)
			ofs_Sliht << log10(Re) << '\t' << 10*log10(1000*1.328/sqrt(Re)) << '\n';
		if((Re>Re_l_cr)&&(Re<1e7))
			ofs_Sliht << log10(Re) << '\t' << 10*log10(1000*0.073/exp(0.2 * log(Re))) << '\n';
		if(Re>=1e7)
			ofs_Sliht << log10(Re) << '\t' << 10*log10(1000*0.455/exp(2.58 * log(log10(Re)))) << '\n';
		//CBS.U_inf++;
	};
	ofs.close();
	ofs_Sliht.close();
	return 0;
	try {
		Task T(&argc, &argv);
		//T.LoadTask(taskFileName.c_str());
		//T.LoadTask("ranque.task");
		//T.LoadTask("zhukov.task");
		T.LoadTask("Blasius.task");
		//T.LoadTask("grad.task");
		//T.LoadTask("rocket.task");
		//T.LoadTask("task.config");
		T.Start();
		//PrintPGrad(T);
	} 
	catch (Exception e){
		e.ShowMessage();	
	};
	//char d;
	//std::cin>>d;
	return 0;	
};

